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A method Is presented for performing efficient and stable finite element 
calculations of heat conduction with quadrilaterals using one-point 
quadrature. The stability In space Is obtained by using a stabilization 
matrix which Is orthogonal to all linear fields and Its magnitude Is 
determined by a stabilization parameter. It Is shown that the accuracy Is 
almost Independent of the value of the stabilization parameter over a wide 
range of values; In fact, the values 3, 2 and I for the normalized 
stabilization parameter lead to the 5-polnt, 9-polnt finite difference and 
fully Integrated finite element operators, respectively, for rectangular 
meshes and have identical rates of convergence In the L 2 norm. Eigenvalues 
of the element matrices, which are needed for stability limits, are also 
given. Numerical applications are used to show that the method yields 
accurate solutions with large Increases In efficiency, particularly In 
nonlinear problems. 
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1. INTRODUCTION 


H 3pA*f 

Vf 5 }A*‘p 


«* 


Because of the versatility of finite element methods for treating complex 
geometries and boundary conditions considerable attention has been focused en 
these methods In heat conduction* One drawback of the finite element method 
compared to the finite difference method Is that presently available 
formulations tend to be more time consuming. For example, in comparing 
standard five point or nine point finite difference formulas in two-dimensions 
with the Isoparametric, bilinear quadrilateral with 2x2 quadrature, one finds 
that In nonlinear heat conduction, a substantial amount of time 1$ used to 
perform the 2x2 quadrature within each element so that the latter can be 
markedly slower. 

The purpose of this paper Is to present techniques throu^- which the bi- 
linear, Isoparametric element for two dimensional heat conduction can be used 
with one point quadrature. Special techniques are needed because when single 
point quadrature Is used, the element matrix contains a spurious singular mode 
In addition to the singular mode associated with the constant temperature 
field. For certain boundary conditions, this singular mode leads to singu- 
larity of the assembled system matrix, which prevents It from being Inverted. 
While the singularity is absent In the transient system matrix, the presence 
of the singular modes In the steady-state matrix will lead to oscillatory 
solutions In which nodal temperatures alternate In sign spatially, and the 
growth of this mode can lead to uninterpretable results. This Is true for 
both explicit and Implicit time Integration procedures. 

This singular mode Is analogous to the hourglass modes found in many fi- 
nite difference codes for transient analysis of contlnua [1] and considerable 
efforts have been devoted to the elimination of these modes In both the finite 
difference and finite element literature [2-4]. 
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In this paper a stabilization procedure Is developed for the quadri- 
lateral element for heat conduction based on the techniques developed In [4] 
so that one point quadrature can be used effectively; In effect , the spurious 
singular mode will be eliminated. The procedure Is described for both the 
conductance vector and the element conductance matrices, so that It can be 
used In both steady -state and transient algorithms with explicit and Implicit 
time Integration. As part of this development, the eigenvalues are obtained 
exactly for both Isotropic and anisotropic heat conduction; this should 
facilitate the choice of a maximum stable time step for explicit time , 
Integration and optimal relaxation factors for Implicit time Integration by 

i 

Iterative equation solvers. 

In Section 2, we review the governing equations for linear and nonlinear 
heat conduction along with the finite element approximations as obtained by a 
variational principle, which are similar to [5] except that they are written 
directly for the nonlinear case. Other nonlinear formulations have been given 
In [6] and [7]. The equations for the one-point quadrature, bilinear, 
Isoparametric quadrilateral are given In Section 3 with the stabilization 
procedure. Section 4 compares the finite element equations given here to 
finite difference spatial semi discretizations based on the standard 5-point 
and 9-polnt molecules. An Interesting result Is that when the mesh Is 
regular, these different molecules can be developed by simply varying the 
stabilization parameter. The eigenvalues of the element matrices are given In 
Section 5, whereas the computer Implementation of this one point-quadrature 
for the heat conduction element Is given In Section 6. 

In Section 7 t we present several example problems. The first two 
examples compare the rate of convergence of this element with 1 point 
quadrature and with 2 x 2 quadrature of the quadrilateral to show the minor 
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effect of reduced Integration on convergence. The remaining problems are 
transient and are intended to show the improvements in speed which are 
possible with this element and the difficulties which result when hourglass 
control is not used; both linear and nonlinear results are presented. 
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2. GOVERNING EQUATIONS AND VARIATIONAL (WEAK) FORMS 


We consider a body a enclosed by a surface r with unit normal n which Is 
subdivided Into a prescribed e-surface and a prescribed flux surface r q . 

We use the following nomenclature 


temperature 


• source per unit volume 


heat flux 


density 


specific heat 


convective heat transfer coefficient law 


■ linear conductivity matrix (k^ • kd^ for Isotropic heat 
conduction) 

The governing equations are: 


- q.j S - PC e 


In Q 


0 « 0 


on r« 


-q 1 n 1 + h(8) - q 


on r 


8 ■ 0, 


In n when t ■ 0 


Standard Indicia! notation is used with repeated subscripts Implying a 
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summation. Hera a comma designates a partial derivative with respect to , 
and a superposed dot designates the time (t) derivative. 

The completion of Eqs. (1) to (4) also requires a heat-i ! aw 


4 


* 


* 


# 


a. 


■ Tj(e» ®»^) (5) 

which for linear heat conduction can be written as 

■ -K,j «.j («) 

The variational or weak form of Eqs. (1) to (4) as given In [8] Is 

fo(«» v) ♦ r(e. v) ■ f(q» s, v) (7) 


where v Is the test function and 


m (e « v) ■ / pc e v do 

a 


•r(0„ v> ■ V M q. da 

a 

f(q, s, v) * / [ q • h(e)]vdr ♦ / s v da 

r a 

q 


( 8 ) 

( 9 ) 

do) 


The finite element equations are obtained by approximating the test func- 
tions and the approximate solution for e(x, t) (trial functions) by shape 
functions Nj. These shape functions are defined In each element, and the 
approximation In each element Is given by using a local separation of 
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where 6? ara the nodal values of the temperature and NODELE Is the number of 
nodes In the element. An Identical expansion Is used for the test function 
v(x) and the space discretization Is performed separately. 

The finite element semi discretization yields the following system of 
ordinary differential equations for heat conduction 

N J + r • f (12a) 

«(0) ■ ^ !12B) 

where M Is obtained from the element matrices by the standard matrix 
assembly of finite elements and r and f are obtained from the element matrices 
by vector assembly. The element matrices are given by 



( 15 ) 
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The above equations are* applicable to both linear and nonlinear heat 
conduction, for r and f may be nonlinear In 9. When the conductance Is 
constant, Eq. (12a) can be replaced by 


M 5 + K e ■ f 


(16) 


where K Is here called the global conductance matrix, which Is assembled from 
element conductance matrices given by 

£ • CKjjf • ! N It1 N J(J <a (17) 

sr 


8 
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3. QUAOR I LATERAL WITH ONE-POINT QUADRATURE OF P° oR 

<* 

TTf? ships functions for a quadrilateral element are written In a 
reference plane $,n In the form 

Hj .^(1 + ?, 0(1 +nj n) (18) 

• 

where Cj, n j are the 5 , n coordinates of node I. If one point quadrature Is 
used, the Integrals In Eqs. (13-1A) and (17) can be computed by simply 
evaluating the Integrands at c ■ 0, n ■ 0 and multiplying by the area. I.e. 
for any function, one-point quadrature gives 

/ f(s, n) da - A f (0* 0) (19) 

a E 

* 

where A Is the area of element E, 

t 

The following equations then hold on the element level 

S”*M E < 20 *> 

* 

r E (D . iVo.O) (20#) 

and the associated element conductance matrix for linear heat conduction Is 

K E U) - l b T D B (21) 

A 

where the superscript 1 designates one-point quadrature. Here 
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(Z2«) 


(22t»> 


(22c) 


The a raa of the element, A, Is given by 


A ’ 7 (*31 y 4Z * *24 y 3l' 


(23a) 


and th* vectors b, are given by 


Si ■ 7 Cy 2 4 *31 *42 *13^ 


Jjj ‘ 7 ( X 42 *13 *24 *31^ 


*IJ ■ *1 - *J 


yw ■ y i - 


(235) 
(23c ) 
(23d) 


For the purpose of Identifying the spurious singular mode of K and Its 
control, we will define two additional column vectors 
T _ r, 1 1 n 


(24*) 


*» iW M 


ORIGIN*- Wj * * 
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h T » Cl. -1. 1. -1] 


(24b) 


and note that 



. o 


bj h ■ 

*0 J +* 


s T h 


(24c) 


These vectors, b^, s and h, span the 4 dimensional space of element nodal 
temperatures and ere shown for a typical quadrilateral In Fig. 1. 

The linear relationship between nodal sources r and nodal temperatures s 
for an element can be written as 

r E • /(D 5 e (25.) • 

■S (b, bj ) S E (255) 

$ 

If we let « E ■ s or e E ■ h, , the orthogonality properties, Eq. (24c ) 

.m 

Immediately lead to the result that r 6 * 0 . Therefore, these two sets of 
nodal temperatures correspond to singular modes of the element matrix 
JC E (*). The first, e E « s , Is expected and necessary since It corresponds to 
a constant temperature field; If a stiffness does not give r E ■ 0 for this 
mode It will not be convergent. We will call this the proper null -space 
of K E W. The second, e E » h, is undesirable and often Is called a spurious 
singular mode, since It can lead to singularity of the assembled finite 
element equations. The presence of an additional singular mode Is often 
called a "rank deficiency" of the element matrix. Note that the two vectors 
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h and s span the null-space of tha element matrix. 

To allmlnata this singular moda, we augment tha element conductance 
matrix by a stabilization matrix [9]. 

JS E * * E(1) + it.b < 2 «) 

where the stabilization matrix Is given by 

.Sitae ■ « 3 I T ' (2?) 

* 

Tha cholca of tho constant l will ba described later. 

This stabilization matrix Is obtained by defining an additional 
generalized thermal gradient g and flux q by 

9 - X T e E (28) 

q ■ « g (29) 


This generalized gradient and flux are added to compensate for the 
contribution to r (8, v) which Is lost due to one point quadrature, so In 
effect we now have Instead of £q. (9) that on an element level 


Thus the element nodal sources are given by 






( 32 ) 


and Eq. (27) ’follows Immediately from (32) and (28). 

The form of x will be chosen so that the following conditions are met: 

I. for any vector of nodal displacements which is defined by a linear 
(or constant) temperature fleld.g - 0 In Eq. (28); 

II. for any other set of nodal temperatures , g # 0 . 

To put this into more precise terms, we designate the vector space of 
nodal temperatures of an element by R* and the null -space of x by rJ*. Since 
the 4 vectors bj, bg, s and h are linearly Independent, they span R*. 

As x I* In R* , we can expand It In terms of these base vectors as follows 

X - *l&l va 2£2 + a 3* + a 4i! (33) 

An arbitrary linear temperature field Is given by 

« (x, y) ■ c x x ♦ c 2 y + c 3 (34) 

and substituting In the nodal values we obtain the following expression for 
nodal temperatures 

S E ■ <1* + «2* + «3 i < 3S4 > 
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1 ,V 

.. '7 

•"' f r 

v { ■ 5 


X T ■ (x lt x 2 , X 3 , X 4 ) 

(35b) 

xj • (y^» >2» y 3 . >4) 

(35c) 

We note the following Identity from [4] 


xj bj - A 1 x 2 a y) 

06 ) 

which Is easily verified by simply substituting In the values of and A. 

The first condition then requires that g (given by Eq. (28)) must vanish 
for all l.e. 

*0 

(a x bj + a 2 bj + a 3 s*+ a 4 ]J) (c x x + c 2 + c 3 s) *0 

for all c.j 

(37) 

Using the orthogonality of s and h and, their orthogonality with b^, and Eq. 
(36) then yields 

1 * J CA £ - (h T x) bj - (h T j£) b 2 ] 

(38) 


Ue will call the vector s^ the proper null -space of R*; Its complement Is of 
dimension 3. 

Since i Is linearly Independent of b^, the 3 together must span the entire 
complement of the proper null -space of R*, so the second condition Is 
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satisfied. 


It Is of Interest to note that the complements of the null -spaces 
of and of J<[ tab (the litter coincides with that of ^ ) are not 
exclusive; l.e. the Intersections of those spaces Is -not empty. This means 


that K 


£stab 


will affect the solution if it Is not linear and the elements are 


not rectangular. Neverthelef j, the stabilization matrix does not affect 
linear or constant fields, so It should not deleteriously affect convergence; 
though this remains to be proven, th»? numerical results In Section 7 confirm 
this fact. 

* 

The stiffness matrix with the stabilization can be written as 


4>., 


L b! ♦ c h h T 
1j - J ~ ~ 


where 


*1j * k 1j + j ) 


k ii<** 1 ❖ 


(41a) 


k < 4 + # 


where t and i are the lengths of the sides of the element. 

A j 


(41B) 


15 


4. 


t- t 3*# # ’ iW s/ r % 

mm<g * 
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COMPARISON WITH FINITE DIFFERENCE FORMULAS AND OTHER FINITE ELEMENTS 


For rectangular and square arrangements of meshes, It Is possible to 
compare this finite element with standard finite difference formulas and fully 
Integrated quadrilateral finite elements. These comparison help in assessing 
the role of the stabilization parameter c and the lack of sensitivity of 
solutions to Its value. 

Fbr a square finite difference mesh, the 5 point and a 9 point mdlecule 
£10] are given In Table 1. The complete stiffness (1 -point quadrature plus 
stabilization) Is also given for c • 3 and * ■ 2. It can be shown by the 
simple assembly of the finite element equations that 

I. c • 3 corresponds to the 5-polnt molecule 

II. c • 2 corresponds to the 9-polnt molecule 

It can also be shown that as defined In Eq. (41), the value of c * 1 
gives the fully Integrated finite element stiffness; while c * 0 of course 
corresponds to the 1 point quadrature stiffness. Thus commonly used finite 
difference and element formulas are associated with a large range of c values. 
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5. EIGENVALUE ANALYSIS OF ELEMENT 

W e consider tht following form of the eigenvalue problem 




k e e E • x E M 2 e E 

tm *0 rm m 


( 42 ) 


where M 6 Is the lumped element capacitance matrix, which Is given by 


M E ■- pc I 

m w 


(43) 


where £ Is the Identity matrix. The system Is associated with the eigenvalue 
problem 

K % ■ \ He (44) 


and according to [12], the largest eigenvalue of any Individual element will 
bound the maximum frequency from above, so 


* max < 


max 

for all E 


*max 


(45) 


Since the stability of Euler Integration requires that 


at < 


max 


(46) 


a time step chosen by 


at 


min 

for all 
E 



(47) 


& 
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In ordeh to obtain the eigenvalues for Eq. (42), we note that the element 
stiffness as given by Eq. (39) Is the sum of 2 terms and the eigenvectors of 
the two terms can be constructed as follows: 

1. The two nonzero eigenvectors of the first term on the right hand side 
of Eq. (39) are linear combinations of bj and b 2 , and since bjand b 2 are In 
the null -space of the second term, this will also be an eigenvector of K E . 

2. The nonzero eigenvector of the second term Is h and, since h Is In 
the null -space of the first term, it Is an eigenvector of K Z . 

The maximum eigenvalue of Eq. (42) can be then be shown to be given by 

4, ■ — (X + Y t /((X- Y) 2 + 4Z 2 ) , 16 c A 2 / k} (48) 

A 

where 


a ■ k/pc 

X * *1J B 11 B lj 

Y * R 1j B 21 B 2j 

Z " **1J B 11 B 2j 


(4$) 


The following special cases are of Interest: 

1. If the material Is Isotropic and the element rectangular 



HIMA | Wm 

*m1n 

where t m ^ n Is the minimum element length, provided that 


, # 
9 



(51) 


where r Is the ratio of the lengths of long side to the short side. This, 
with Eq.(47), yields the following condition for stability 


at < 



(52) 


3 

11. For square meshes with a distance i between nodes and with c > - , the 

2 

maximum eigenvalue is given by the second term In Eq. (48), 1. c. 



(53) 


Remark 1. The eigenvalue In Eq. (53) governs the time step for the 5-polnt 
and 9-polnt difference formulas (c ■ 3 and 2 respectively), so the stable 
time step for these difference formulas Is smaller than for the finite element 
method. This contrasts with the findings In Cll] and [12], where the opposite 
was found because (1) less accurate bounds were used for the elgenvaules and 
(2) the consistent capacitance matrix was used. 

.. .Remark 2. The stability limit for the time step resulting from Eq. (53) for 
the 5-polnt difference formula (« ■ 3), agrees exactly with the result of a 
Neumann analysis given In (13) 
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and Hourglass Control 

For Quadrilaterals 

For simplicity, we have dropped the superscript E In this section. We 
first define explicitly the one point quadrature element; vector r^ 1 ^ and the 
stabilization element vector r h employing hourglass control. Then the one 
point quadrature with hourglass control element vector r Is equal to the sum 
of r and r^. 

One Point Integration 

As give In eq. (14), the element vector £ Is: 

• 

r I ■ * la <»l.x V N l,y V « 

(55) 

where q., and q are (see Eq. (22c) ): 
* y 


q x ■ - (k H 9 x + k 12 9 y> 

(56a) 

q y * - ( k 12 9 x * k 22 9 y> 

(560) 

and 

* 


9 x * C y 24 * e l " 0 3> * y 31 (e 2 * ®4 )3 

(see) 

9y * C x 42 (0 1 * ®3 } * X 13 ( ®2 * ®4 )] 

(566) 

Employing one point Integration 
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r l 1> * 5 ll q x + 6 2I q y * “tl q 1 


(57) 
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where and b 2I are defined In eas. (23b) and (23c) respectively. 
Hourglass Control 

The stabilization vector as defined in eq. (32) is: 



where 

m m 

q ■ c g 


(58) 



(59) 


1 • 7 C («! - * « 3 * »4> ' <V* + Vy )] 

a x * X 1 “ X 2 * x 3 ' *4 
a y * *1 " y 2 * y 3 * y 4 

The ncdal components of this hourglass vector are: 


(60a) 

(60b) 

(60c) 


p l" * h I^ * l b 1l x iJ h J^ 

Therefore the element vector using one point quadrature and hourglass 


control for a quadrilateral Is: 


r[ .h, q- I D,,*, jhjq + b,^ 
After some algebra one can show that 

P I • b XI< + °2I^y + h I< 


( 62 ) 


. ( 63 ) 


where 


7. Numerical Results 


* 



A two dimensional finite element pilot computer code Incorporating the 
methodologies described In the previous sections has been written to evaluate 
the performance of this one point quadrature element and our critical time 
step estimates. Four numerical examples are presented to demonstrate the 
accuracy, stability criterion and efficiency of these proposed methods. 

e 

Results are compared with exact solutions or approximate solutions using two 
by two quadrature. All computations are performed on a COC Cyber 170/730 
computer In single precision (60 bits per floating point word). For the 
transient analysis, a lumped capacitance matrix Is used and the predictor- 
corrector explicit algorithms with a ■ 0.5 used In [8] are employed to carry 
out the time Integration. 

Example 1: Convergence Study of a Unit Square Plate with Prescribed 
Temperatures 

Due to symmetry of the geometry and prescribed temperatures, only half of 
the unit square plate Is modelled with 32 (4 x 8), 128 (8 x 16) and 200 (10 x 
20) elements respectively. These three finite element meshes are depicted If? 
Fig 2. Side BC Is a line of symmetry (Insulated). Sides AO and DC are 
prescribed a constant uniform temperaure of 0.0, while sides AB Is prescribed 
with a constant temperature distribution of sin nx where the x-axIs Is defined 
by Joining node A to node 8. Hence the temperatures at nodes A and B are 0.0 
and 1.0 respectively. The exact steady-state solution Is given by: 

e exact m s1nh n (i a Q«y ) sin nx/slnh it 

Two values of the stabilization parameter c were tested. For c ■ 1.0, this 
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element Is Identical to the two by two quadrature element (since the elements 
are rectangular)} whereas for c • 0,0, It Is Identical to the one point 
quadrature element. The temperature profiles along 8C obtained from these 
three finite element meshes (with c « i.0) and from the analytical solution 
are also desplcted In Fig. 2. The finite element solutions of the case 
c • 0.0 differ from those of the case c ■ 1.0 In the third or fourth 
digit. Therefore they are not plotted. As can be seen, the finite element 
solutions are virtually Identical to the exact solution. 

We also computed the L 2 error norm for these solutions as follows: 

« 

a 

« 

E . ( / A • V* 

J 

where t ■ e #x * ct - and A Is the area. The total L 2 error, E, Is computed 
using a 5 x 5 quadrature In each aliment. We uotalned convergence rates of 
1.899, 1.908 and 1.930 for the cases of c * 2.0, 1.0 and 0.0, respectively, 
which agree reasonably with the theoretical convergence rate of 2. 

Remark 1. The reduction of the quadrature rule from 2 x 2 to 1 has no 
significant effect on the convergence rate. 

Remark 2. It Is possible to solve this problem with the stabilization 
parameter c » 0 because the boundary conditions eliminate the rank deficiency 
of the assembled mesh. This Is not always possible, as will be seen 
subsequently. 

Remark 3. The convergence rate of the 9 point Laplaclan (c ■ 2.0), which has 
a much smaller truncation error, shows no Improvement over the finite element 
method . 

Example 2 : Convergence Study of a Circular Plate with a Heat Source 

Due to double symmetry, only a quarter of the circular plate (which Is 
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heated with a uniform constant heat sourca ■ 1.0} is modelled with 12, 48 

* 

and 192 elements respectively, Tht finita elements mashas ara shown In Fig. 

3. It should ba observed that soma of the quadrilateral elements ara quite 
skewed. The exact solution for this circular plate with radius r ■ 5.0, 
thermal conductivity k » 0.04 and a constant temperature of 0.0 at r ■ 5.0 Is 
given as: 


0 Exact ^ „ 6f25 ( 25 . p 2) 

As In the preceding example, two values of the stabilization parameter of 1.0 
and G.O respectively are tested. However, due to the skewness of the 
elements, the c * 1.0 elements are not the same as the two by two quadrature 
elements. The c ■ 0.0 elements are still Identical to the one point 
quadrature elements. The temperature profiles along nodes 1 to 45 obtained 
from these three finite element meshes (with t ■ 1.0) and the exact solution 
are also desplcted In Fig. 3. Again, we found that the finite element 
solutions of the case c * 0.0 differ from those obtained using c * 1.0 In the 
third or fourth digit. Therefore they are not plotted. The polntwlse 
convergence of this stabilized element Is cleared shown In the plot. 

We obtained convergence rates of 1.955 and h924 for the cases 
of c • 1.0 and 0.0 respectively which agree well with the expected convergence 
rate of 2.0. 

Example 3 : Linear Transient Thermal Ar.alysls of a Wedge 

The problem statement Is depicted on the top of Fig. 4. The finite 
element mesh consists of 100 elements and 121 nodes. The thermal dlffusivlty 
of the wedge Is 0.001. The Initial temperature for all the nodes Is 0.1. All 
four sides are Insulated. The heat load which Is also shown in Fig. 4 is 


-ss&ss* 

applied at node l. A constant time step of 1.0 Is used for this problem, 
fhls time step Is computed according to Eq. (48)* The temperature-time 
histories at four different locations are presented also in F1g\* 4. These 
results are obtained using c * 1.0 stabilized element. These results are 
virtually identical to those obtained using two by two quadrature elements* 

For values of c ■ 0.8 and 1.2. the peaks at node 1 are about 4% below and 4% 
above the solution with 2x2 quadrature, therefore c * 1.0 Is recommended. 
Example 4 ; linear and Nonlinear Transient Thermal Analysis of a 
Circular Plate 

The N med1um N finite element mesh (48 elements with no heat source) shown 
In F1g.^1s employed for this problem. The heat load which 1 $ shown In Fig. 4 
Is applied at node 1. The Initial temperature for all the nodes Is 0.1. All 
boundaries ere Insulated, The thermal dlffuslvlty of the plate Is 0.004. 
According to Eq. 48, It corresponds to a critical time step of 1.0. Two 
hundred time steps are run (at the critical time step) to obtain the 
temperature-time histories shown In Fig. Sa. These results are obtained 
using c • 1.0 and the solutions are virtually Identical to those obtained 
using twu by two quadrature elements. However, we obtefned severe spatial 
oscillatory solutions using c ■ 0.0 ror this problem (see Fig. Sb). 

In order to demonstrate the effectiveness of this one point quadrature 
element with stabilization, the thermal dlffuslvlty, a, Is changed to: 

a ■ 0.004 (1.0 + O.Ole) 

to make the problem nonlinear. A constant time step of 0.2 Is used for this 
nonlinear problem. The computed solution using c ■ 1.0 are also presented In 
Fig. 5a. These results are almost the same as those using two-b,y-two- 






quadrature elements except there is a 3* difference in the peak temperature of 
node 1. However, we gain a factor of 4.38 In solution time by employing the 
stabilized one-point element as compared to 2 x 2 quadrature. Although a 
factor of 4.8 would be expected, the savings are actually greater because the 
shape functions need not be evaluated at quadrature points In this procedure. 

t 

8. CONCLUSIONS 

In this paper, an efficient computational method has been developed for 
the linear and nonlinear heat conduction with a quadrilateral element. A 
computationally -useful method of estimating the critical time step for this 
element in explicit time Integration Is given. The computer implementation 
aspects as well as the evaluation of the performance of this new element as 
applied to two-dimensional steady and transient thermal analysis are also 
presented. 

Numerical results show: 

(1) this method yields accurate solutions, 

(2) the great Increase In computational efficiency especially In nonlinear 
analysis, and 

(3) the Importance of this method as applied to three dimensional and/or 
nonlinear thermal analysis. 

Comparison with flnlta difference formulas has shown that various values 
of the stabilization parameter, the 5-polnt and 9-polnt molecules can be 
obtained. The convergence rate, however, appears to be Independent of c which 
means It Is Independent of the order of quadrature In the finite element 
method. 
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TABLE 1 


Relationship between stabilization parameter e and finite difference formulas 
square mesh 
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9 point, Ref. [10] 
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FIGURE CAPTIONS 


Fig. 1 Base Vectors for a typical quadrilateral. 

Fig. 2 Convergence study of a Unit Square plate with Prescribed Temperatures . 
Fig. 3 Convergence study of a circular plate with a heat source. 

Fig. 4 Linear transient thermal analysis of a wedge. 

Fig. 5 Linear and Nonlinear Transient Thermal Analysis of a Circular 
Plate (a) results obtained with the stabilization matrix 
(b) oscillatory solutions obtained without the stabilization matrix. 




















